Preprint typeset in JHEP style - PAPER VERSION 



NLO Higgs boson production via gluon fusion 
matched with shower in POWHEG 



Simone Alioli 

Universita di Milano-Bicocca and INFN, Sezione di Milano-Bicocca 
Piazza della Scienza 3, 20126 Milan, Italy 
E-mail: Simone . AlioliOmib . inf n . it 

Paolo Nason 

INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milan, Italy 
E-mail: Paolo.NasonOmib.infn.it 

Carlo Oleari 

Universita di Milano-Bicocca and INFN, Sezione di Milano-Bicocca 
Piazza della Scienza 3, 20126 Milan, Italy 
E-mail: Carlo . Dlear iOmib . inf n . it 

Emanuele Re 

Universita di Milano-Bicocca and INFN, Sezione di Milano-Bicocca 
Piazza della Scienza 3, 20126 Milan, Italy 
E-mail: Emanuele.Re@mib.infn.it 

Abstract: We present a next-to- leading order calculation of Higgs boson production via 
gluon fusion interfaced to shower Monte Carlo programs, implemented according to the 
POWHEG method. A detailed comparison with MCONLO and PYTHIA is carried out for several 
observables, for the Tevatron and LHC colliders. Comparisons with next-to-next-to-leading 
order results and with resummed ones are also presented. 



Keywords: QCD, Monte Carlo, NLO Computations, Resummation, Colhder Physics 



Contents 



1. Introduction 1 

2. Description of the calculation 2 

2.1 Kinematics 2 

2.1.1 Born kinematics 2 

2.1.2 Real-emission kinematics 3 

2.1.3 Inverse construction 4 

2.2 Cross sections 5 

2.2.1 Born contribution 5 

2.2.2 Virtual corrections 5 

2.2.3 Real corrections 5 

2.2.4 Collinear remnants 6 

3. POWHEG implementation 7 

3.1 Generation of the Born variables 7 

3.2 Generation of the radiation variables 11 

4. Results 12 

4.1 POWHEG - MCONLO comparison 12 

4.1.1 Tevatron results 13 

4.1.2 LHC results 14 

4.2 POWHEG - PYTHIA comparison 15 

4.3 The Pt distribution in POWHEG 19 

4.4 Next-to-leading logarithmic resummation 24 

5. Conclusions 26 



1. Introduction 

Gluon fusion is the dominant Higgs boson production mechanism both at the Tevatron 
and at the LHC. Radiative corrections to this process are known to be large [1, 2, 3], and 
it is thus important that shower generators that do include them are made available to the 
experimental collaborations. In fact, one such generator already exists, namely the MCSNLO 
implementation [4] of Higgs boson production. 

In this work we present a next-to-leading order (NLO) calculation of Higgs boson 
production via gluon fusion, interfaced to shower Monte Carlo programs according to the 
POWHEG method. Unlike the MCSNLO implementation, our generator produces events with 
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positive (constant) weight, and, furthermore, is not tied to the HERWIG shower Monte Carlo 
program. It can be easily interfaced to any modern shower generator and, in fact, we show 
results of POWHEG interfaced to HERWIG [5, 6] and to PYTHIA [7]. 

The POWHEG method was first suggested in ref. [8]. In ref. [9] a detailed general de- 
scription of its application to collider processes was given. Until now, the POWHEG method 
has been applied to ZZ pair hadroproduction [10], heavy-flavour production [11], e^e~ 
annihilation into hadrons [12] and into top pairs [13], and Drell-Yan vector boson produc- 
tion [14, 15]. We have built our implementation of the Higgs boson production by following 
closely the formulae and results of ref. [9] . 

Much of our phenomenological section will be devoted to study the comparison of our 
result with that of MCQNLO. We find fair agreement between MCONLO and POWHEG results, 
except for the p-r distribution of the Higgs boson, and consequently of the hardest jet, 
in the high-p^ region. In this region, the POWHEG distributions are generally harder. We 
have shown that this is due to next-to-next-to-leading order (NNLO) effects in the POWHEG 
formula for the differential cross section. We checked that these effects actually bring our 
result closer to the NNLO one [16]. Other relevant discrepancies are found in the rapidity 
difference of the Higgs boson and the hardest jet. The dip produced by the MCONLO program, 
found in previous implementations [10, 11, 14], is present also here. We remark that this 
seems to be a general feature of MC@NLO, since other calculations do not find effects of this 
kind [17, 18, 19]. 

The paper is organized as follows. In sec. 2 we describe how we performed the calcula- 
tion for the Higgs boson cross section at the next-to-leading order. In sec. 3 we discuss the 
POWHEG implementation. In sec. 4 we show our results for several kinematic variables and 
compare them with the MCSNLO [4] and PYTHIA 6.4 [7] shower Monte Carlo programs. A 
comparison with next-to-next-to-leading order results, as well as with analytical resummed 
ones is also carried out. In sec. 5, we give our conclusions. 

2. Description of the calculation 

In this section we fix our kinematic notation, and give the Higgs boson production differ- 
ential cross sections up to next-to-leading order in the strong coupling Og. 

2.1 Kinematics 

2.1.1 Born kinematics 

The Born process has a single partonic contribution, gg — ^ H. Following the notation of 
ref. [9], we denote with fc^ and fcg the incoming gluon momenta, aligned along the plus 
and minus direction of the z axis, and by k\ the outgoing Higgs boson momentum. If Kf^ 
and Kq are the momenta of the incoming hadrons, then we have 




(2.1) 



where the momentum fractions, and momentum conservation reads 



(2.2) 
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We introduce the Higgs boson invariant mass squared and rapidity 



M' = kl ^ = ^l°gKI' (2-3) 

1 1 



so that the set of variables *i = |M^,y} fully parametrizes the Born kinematics. Prom 
them, we can reconstruct the momentum fractions 



/^^ Y -Y 

Xe = y -g-e , = Y , (2.4) 

where S = {K^ + Kq)'^ is the squared center-of-mass energy of the hadronic collider. The 
Born phase space, in terms of these variables, can be written as 

We generate the Higgs boson virtuality according to a Breit-Wigner distribution, i.e. we 
make the replacement^ 

d[M —mij) o • (2.6) 

The decay of the Higgs boson is left to the shower Monte Carlo program, since, being the 
Higgs boson a scalar, no spin correlation can arise. 

2.1.2 Real-emission kinematics 

The real emission processes have an additional final-state parton, so that momentum con- 
servation reads 

fee + fee = fci + A;2 , (2.7) 

where ki is the Higgs boson momentum and k2 is the momentum of the additional final- 
state parton in the laboratory frame and 

k^ = x^K@ . (2.8) 

Since we regularize the infrared divergences in the Prixione, Kunszt and Signer (FKS) 
subtraction scheme [20, 21], we introduce the appropriate set of radiation variables. In the 
partonic center-of-mass frame, the final-state parton has momentum 

^2 = ^2° (IjSin^sin^, sin0cos(^,cos0), (2.9) 

and we use the set $rad = {^^U^'P} as radiation variables, where 

A^2° = ^C, y = cos0, (2.10) 



^In order to compare our result with other programs, we have also used slightly different forms of the 
Breit-Wigner distribution, that will be illustrated in due time. 
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and 

s = {k^ + k^f = ^ (2.11) 

is the partonic center-of-mass energy squared. Since there are no final-state coloured par- 
tons at the Born level, we have to deal with initial-state singularities only. The soft 
singularity is characterized by ^ ^ 0, while the collinear limits (/c2 parallel to the ® or e 
incoming directions) are characterized by y ^ 1 and y — > — 1 respectively. 

2.1.3 Inverse construction 

The set of variables $2 = {M'^,Y,^,y,(l)^ fully specifies the real-emission kinematics. In 
fact, given these variables, we can reconstruct all the momenta. Using eq. (2.4), we can 
compute the underlying Born momentum fractions x@ and, following sec. 5 of ref [9], we 
have 

" " V 2 - ^1 + y) ' " " VT^\I 2 - ^1 - y) ' ^'-'"^ 

with the kinematics constraints 

< c < Uy) , (2-13) 

where 



^M(y) = 1 - max 



2(l + y)a:^ 



V(i + xiy{i - y? + iQyxl + (1 - y)(l - xl) ' 

2(1 -y) 4 ^ 

V(l + xlY{l + yY - 16y x| + (1 + y)(l - xl) 



(2.14) 



The momentum of the final-state parton in the partonic ccnter-of-mass frame is given by 
eqs. (2.9) and (2.10). We then make a longitudinal boost from the center-of-mass frame 
back to the laboratory frame, with boost velocity 

/3=^^^, (2.15) 

Xq + Xq 

to obtain k2 from ^2 

k2=^Lk2. (2.16) 
Prom momentum conservation, we reconstruct the Higgs boson momentum 

ki = x^K^ + XqKq - k2 . (2.17) 

Finally, the two-body phase space can be written in a factorized form in terms of the Born 
and radiation phase space 

= dx^ da;e(27r)^5^(A;e + k^ - ki - ^2) (^27^)32^? (27r)32fc^ " ^^"^ ' ^^'"^^^ 

where 

d^rad = ^ dy d<t> . (2.19) 
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2.2 Cross sections 



In order to apply the POWHEG method, wc need the Born, real and virtual contributions 
to the differential cross section, i.e. the squared amplitudes, averaged over colours and 
helicities of the incoming partons, and multiplied by the appropriate flux factor. 

2.2.1 Born contribution 

At Born level, Higgs boson production via gluon fusion proceeds through the coupling of 

the Higgs boson to a heavy-quark loop. The squared matrix clement for the lowest-order 
contribution, averaged over colours and helicities of the incoming gluons, and multiplied 
by the flux factor 1 / (2M^ ) , is given by 



B. 



99 



576^2 



^^rQ[l + {l-TQ)f{TQ) 



(2.20) 



where tq = Am'q/M'^, and the sum runs over the heavy flavours with mass mq circulating 
in the loop. The function / is given by 



firq) = <^ 



arcsm 



log 



1 + V(l-rQ) 



rQ>l, 
rq < 1- 



(2.21) 



1- va-TQ)^ 

In our implementation we only retain the contribution coming from the top quark. 
2.2.2 Virtual corrections 

In the calculation of all NLO corrections, we have used an effective Lagrangian, where the 
heavy-quark degrees of freedom have been integrated out. This corresponds to take the 
mq CO limit. 

We have regularized the infrared divergences according to the conventional dimensional 
regularization method, i.e. we have set the space-time dimensions D = A — 2e. 

The finite soft-virtual term, obtained from the sum of the divergent virtual contribu- 
tions and of the integral over the radiation variables of the counter-terms is given by (see 
eq. (2.99) of ref. [9]) 



~ 27r 



11 
3" 



^C.-^Tpn/) logg + :^C. + ^C. 



11 

y 



2_7r2 
3 



B. 



99 ■ 



(2.22) 



In deriving this equation we have set = 1- We indicate with /j,^ and Hp the renormaliza- 
tion and factorization scales, respectively. 

2.2.3 Real corrections 

At NLO, there are four subprocesses that contribute to Higgs boson production: gg Hg, 
gq Hq, qg Hq and qq Hg, where q runs over all possible quark and antiquark 
flavours and q and q are conjugate in flavour. The respective squared amplitudes, averaged 
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over the incoming helicities and colours and multiplied by the flux factor l/(2s) are given 

by 

_ al Gp 1 s^ + t^ + u^ + M^ 

~ Utt V2 2s stu ' ^ 

^ al Gp 1 s^ + v?' /„„,x 

'=« = -2^7f2;— • 

,e,,._^^lf!±^, (2.25) 

277r ^ 2s « ' ^ ' 

T^m = ^-k — — — ' (2-26) 



Slvr V2 2s s 



where 



s = (A;e + A;er = t = - ^2)' = - ^ C (1 + y), ^ = (^® - ^2)' = "2 C (1 - y)- 

(2.27) 

In terms of the FKS variables we then have 

= tJ ^ I + + + 1' + - ?(T^ • 

where the singular behavior for a soft — 0) or collinear gluon {y — ±1) is clearly 
manifest. Notice that the contribution Tlqq is not singular and has no underlying Born. 

2.2.4 Collinear remnants 

After the subtraction of the initial-state collinear singularities into the parton distribution 
functions, finite collinear remnants are left over. The kinematics of these terms is Born- 
like. More precisely, we can introduce two sets of variables, = {M^,y, z}, such that 
momentum conservation reads 

zx^K^ + XQKQ = ki (2.32) 

for the © direction and 

Xe-fC® + zXqKq = ki (2.33) 
for the e one. We can then associate an underlying Born configuration $1 such that 

= zx^Kq , /cq = XqKq , ki = ki (2.34) 

for the © direction, and 

^® ~ ^e-^e ) ~ ^^e-^e ) ki = ki (2.35) 
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for the e one. 

The collinear remnants arc given in eq. (2.102) of rcf. [9], where we have fixed = 1 
and 5i = 2 and chosen the MS rcnormalization scheme. For the © direction and for the 
two different real-term contributions, they are given by 



ffS'(*.,e)=£cJ(l-.)i±(i^ 



1 



1-z 
+ 2 



log 



log(l - z) 



1-z 



+■ 



z + - ^- + z{l-zf 



l-z 



+ 2 



log 

log(l - z) 
1 - z 



+ Z^Bgg, (2.36) 
M2 



Bgg. (2.37) 



The other two collinear remnants, ^0*^(^1,0) Qq' {^1,0)7 have the same functional form 
of Q^^ Q^^ (^1,®) respectively, since Bgg only depends upon kf. 



3. POWHEG implementation 

3.1 Generation of the Born variables 

The first step in the POWHEG implementation is the generation of the Born kinematics. 
According to ref. [9], we introduce the -B(^i) function, defined as 

= Bgg{^l)+Vgg{^i) + ^ d^.^A^^Rgg {^1,^,^) 

+ Yl [^19 ^rad) + Rgq ^rad) 

9 



GSn*i,e) + E^I'(*i.e) , 

(3.1) 



dz 

+ I — 



G|^(*i,e)+E^e(*i-e) 
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where 



Gr(*i,e 
G|^(*i,e 



Sggi^l) Cgg{x^,Xe), 

Vgg{^l) CggiXg^Xe), 
Ugg{^l,^rad) i:gg{x<s,XQ) , 

Ugq{^l,^r3d) ^gq{x@,XQ) , 
'R-qg{^l,^rad) ^qg{XQ,XQ) , 



99 



Ql'{^l,^)Cqg{^,X 



= a|^(*i,e)>C,p(xe,^), 



91 



(3.2) 
(3.3) 
(3.4) 
(3.5) 
(3.6) 

(3.7) 
(3.8) 
(3.9) 
(3.10) 



with X(^,Xq given in eq. (2.12) and the luminosity C is defined in terms of the parton 
distribution functions /®(x@,/Xp) 

Cffi{x^,XQ) = ff{x^,i4) ff{xQ,ixl). (3.11) 

Observe that the TZqq term docs not appear in B, since it does not have a vahd underlying 
Born. It is just generated separately, as described at the end of this section. 

All the integrals appearing in eq. (3.1) are finite. In fact, according the the FKS 
subtraction scheme, the hatted functions 



7^. 



1 fl 



ij 



C 12 



1 



1-2/ 



+ 



1 



1 + y 



[(1-2/2) ^2^..] 



(3.12) 



have only integrable divergences. Some care should still be used when dealing with the 
plus distributions. In order to illustrate this, we explicitly show how to deal with the Tlgg 
term, that is the most singular one. According to eq. (3.12), it can be written 



V — ^ 



^99 



127r V2 I 



2 + yl±^l±l^^ + 2(1 - 0' 



1 /I 



2 



1 



1-2/ 



+ 



1 



1 + y 



1 

I 

(3.13) 



Inserting now the expression (2.19) of d^^ad into eq. (3.1), we have 



<^^rad ^gg{x^j Xq) T^gg 



_M2 

(47r)= 



d(f) ^gg{x^, Xq) TZg 



^99 ' 

(3.14) 

where ^M(y) is given in eq. (2.14). The integration over the azimuthal angle (f) is straight- 
forward, giving an overall multiplicative factor of 27r. Considering then the (1/(1 — y))+ 
term only, we get an integral of the form 



dy 



+ 



1 



(3.15) 
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where 



1 



27r 



_ 2 (47r)3 (1 _ 



Recalling the definition of the plus distributions 



/(e)-/(o) 

/(y)-/(i) 



and making the change of variable 



we are left with 



I = 



dy 



-1 



1-y 



fiC^y) ^99 (a;®,a;e) 



1 



+ dy 

-1 1 - y 



/(0,y) logeM(y)-/(0,l) log^M(l) ■^55 (-^e) ^e) 



= f dy [ di—^^\\f{^,y)£'gg{x^,XQ)-f{0,y)£gg{x^,: 

Jo Jo ^ — y^[^ 



1-r 



(3.16) 

(3.17) 
(3.18) 

(3.19) 



+ / z f{0,y) log^Miy) - /(0,1) log^M(l) JCgg{x^,Xe) , (3.20) 

Jo ^ ~ y 

where we have used the expression of x@ of eq. (2.12) and ^m(I) = 1 — (see eq. (2.14)). 
In the last line we have made the further change of variable 



y = -l + 2y, 



(3.21) 



so that all radiation variables are mapped into a cubic unit volume. The integral I is now 
manifestly finite and can be computed numerically. 

The same manipulations should be applied to the z integration of the coUinear rem- 
nants in eq. (3.1). For example, concentrating on the two plus distributions in the Q^^ 
term, we have to deal with integrals of the form 



(3.22) 



dz 



1 - z 



fiz) = \ \og\l - xe)/(l) + dl l^ill^ [/(^) _ /(I)] , (3 
^ Jo 1 - £ 



.23) 
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where f{z) is finite in the z — > 1 hmit and we have made the change of variable 

z = x^+iil-x^) . (3.24) 
At the end of this procedure, the most general form one can obtain for B is 

B{^i) = D{^i) + j\iE[^,,i) + dy j\iF[^i,ly) , (3.25) 

and we can define the function 

=Z}(*i) +F(*i,|,y) , (3.26) 

so that 

B{^i)= j\yj\iB(^^i,lyy (3.27) 

In order to generate the underlying Born kinematics, we first compute the two distinct 
contributions to the total cross section, defined by 



where 



and 



(7h^_ = J d*id$,adi?99-(*l,^rad) , (3.29) 
Rqq{^l,^ra.d) = T^qqi^ I, ^Ta.d) ^qq{Xs, Xq) . (3.30) 



We then decide whether the event is a 5 event or a Rqq one, with a probability equal 
to as/cFtot Bind cth^^ /o"tot respectively. In case of a ^ event, the generation of the Born 
variables #i is performed by using the integrator-unweighter program MINT [22] that, after 
a single integration of the function B{^i,^,y) over the Born and radiation variables, can 
generate a set of values for the variables distributed according to the weight 

B{^i,^,y). We then keep the generated values only, and neglect all the others, which 
corresponds to integrate over them. The event is then further processed, to generate the 
radiation variables, as illustrated in the following section. In case of a Rqq event, one uses 
the same method used for the B case, except that, at the end, one keeps the whole set of 
Born plus radiation variables, that fully defines the kinematics of a real event. In this last 
case, one does not need to do anything else, and the event is passed to the Les Houches 
Interface, to be further showered by the Monte Carlo program. 
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3.2 Generation of the radiation variables 

Radiation kinematics is generated using the POWHEG Sudakov form factor 

A(*i,Pt) = expj- J d^,^^^^^^e{k^-p^)^, (3.31) 
where we have defined 

R ^rad) = Rgg ^rad) + Yl [^19 ^rad) + Rgq ^rad)] , (3-32) 

= (*i) > (3-33) 

and 

kl = le{i-y') = ^^^e{i-y') (3.34) 

is the exact squared transverse momentum of the radiated parton. The factorization and 
renormalization scales in eq. (3.31) should be taken equal to k^, in order to recover the 
correct leading logarithm (LL) Sudakov behavior^. 

To generate the radiation variables, we use the veto method. This requires to find a 
simple upper bound for the integrand in eq. (3.31) 



M2 e i?(*l,$rad) 



(47r)3 (1-^)2 5(*i) 

A suitable upper bounding function is given by 



(3.35) 



where N is determined by spanning randomly the whole phase space and imposing that U 
is larger than the integrand function. The generation of the event according to the bound 
(3.36) is documented in great detail in Appendix D of ref. [23], and we do not repeat it here. 

The POWHEG differential cross section for the generation of the hardest event is given 

by 

da = S(*i) d*i I A(*i,K"") + A{^,,p^) ^ ^l^'^'^j'^^ d<l>rad I 



+ J2 *rad) d*id$,ad , (3.37) 



where the last term in the sum is the non-singular real contribution. In the B and Rqq 
functions, the renormalization and factorization scales, /7,j, and yUp, should be taken of the 
order of the hard scale of the process, i.e. the Higgs boson mass or its transverse mass. 



We will show in sec. 4.4 how it is possible to reach next-to-leading logarithmic accuracy for this particular 

process. 
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During the generation of radiation, the two scales should instead be taken equal to the 
transverse momentum of the produced radiation, in order to recover the correct Sudakov 
form factor. 

We remark that, in the formula for the strong coupling constant used for the generation 

of radiation, we have properly taken into account the heavy-flavour thresholds. That is to 

say, when the renormalization scale /ip crosses a heavy-flavour mass threshold, we change 

the number of active flavours accordingly. Furthermore, as discussed in refs. [9, 10], we use a 

(5) 

rescaled value Amc = 1-569 in the expression for a^, in order to achieve next-to-leading 
logarithmic accuracy in the Sudakov form factor (see sec. 4.4 for more details). 

4. Results 

In this section we present our results, obtained for the Tevatron and the LHC, and the 
comparison done with MCONLO and PYTHIA. We have used the CTEQ6M [24] set for the 
parton distribution functions and the corresponding returned value = 0.226 GeV. In 

the generation of the radiation, we have fixed the lower cutoff of the transverse momentum 
to the value p™"^ = \/5 A^c- The renormalization and factorization scales have been taken 
equal to the Higgs boson transverse mass m!^ = ^ m'jj + [p^Y- 

No acceptance cuts have been applied in any of the following plots. 

4.1 POWHEG - MCONLO comparison 

We have compared our results with MCSNLO, the only existing program where NLO Higgs 
boson production via gluon fusion is merged with a shower Monte Carlo program. Since 
MCQNLO uses only the HERWIG [5, 6] angular-ordered shower, we have also interfaced POWHEG 
with HERWIG, in order to minimize effects due to differences in the shower and hadronization 
algorithms. 

MCSNLO generates the Higgs boson virtuality according to the Breit-Wigner form 

}_ niH Th 

For the purpose of this comparison we have thus used the same form. We have considered 
two different sets of values for the Higgs boson mass and width: mH = 120 GeV with 
Th = 3.605 MeV and mn = 400 GeV with Th = 28.89 GeV. 

Both in POWHEG and in MCSNLO there is the option to retain the full top-mass dependence 
in the Born cross section, i.e. to use a finite tq value in eq. (2.20). We have then the choice 
to generate our Born variables by fixing rrit = 171 GeV in the B term in eq. (3.1) or by 
sending mt — ^ oo. Since we have computed the real-radiation term only in the rrit ^ oo 
limit, we have to use the same limit in the calculation of the Born term in the Sudakov form 
factor (3.31), in order to recover the correct Altarelli-Parisi behavior when the collinear 
limit is approached. 
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(4.1) 
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Figure 1: Comparison between PDWHEG and MCONLD for the rapidity, invariant mass and transverse- 
momentum distributions of a Higgs boson with mn = 120 GeV, at Tevatron pp colhdcr. 




-4 -2 2 4 -4 -2 2 4 



yjet yjet yn 

Figure 2: Comparison between PQWHEG and MCONLO for the rapidity of the leading jet and the ra- 
pidity difference of the Higgs boson and the leading jet, defined according to the SISCQNE algorithm, 
with different jet cuts. 



4.1.1 Tevatron results 

In fig. 1 we show a comparison between POWHEG and MCSNLO for the rapidity, invariant mass 
and transverse- momentum distributions of a Higgs boson with mass mn = 120 GeV, at 
the Tevatron pp collider. The lowest order mt-dependence is retained. A blowup of the 
transverse-momentum distribution near the low-p^ region is also shown. There is good 
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agreement between the two programs, except for the transverse momentum distribution at 
high (we will comment more on this issue in sec. 4.3). 

In fig. 2 we compare the leading jet rapidity and the difference in the rapidity of the 
leading jet and the Higgs boson. The jet is defined using the SISCDNE algorithm [25] as 
implemented in the FAST JET package [26], setting the jet radius R = 0.7 and the overlapping 
fraction / = 0.5. As in previous POWHEG implementations, we notice a dip in the MC@NLO 
jet rapidity distribution, which is enhanced in the difference. We have already extensively 
discussed this fact in sec. 4.3 of ref. [14]. 

In fig. 3, we compare the transverse- momentum distributions of the leading jet, recon- 
structed with the SISCONE and the /cr algorithms (included in FAST JET). A lower 10 GeV 
cut on jet transverse momentum is imposed. The high-p^ discrepancy reflects the same 
behavior found for the Higgs boson transverse-momentum distribution (see sec. 4.3). 



4.1.2 LHC results 

From fig. 4 to 6 we carry out a similar analysis for the LHC pp collider. The difference in 
the hardness of the distributions is more evident here than at the Tevatron. The other 
plots show instead a good agreement between the two codes, apart from the aforementioned 
dip in the leading-jet rapidity distributions. 

We have also made some comparisons with a different value of the Higgs boson mass. 
We have chosen mn = 400 GeV, where the ratio between the Born cross sections evaluated 
with rrit = 171 GeV and — > oo is close to its maximum value and roughly equals 3. The 
results are shown in fig. 7 and 8. We see that, in this case, the dip in the rapidity of the 
hardest jet in MCONLO is extremely marked. 

In the study of ref. [17], carried out in the framework of heavy-flavour production, the 
origin of the rapidity dip was tracked back to an even stronger dip in the pure HERWIG 
distribution, that the MC@NLO correction was not able to properly fill. The same pattern is 
also observed in the present context, as can be seen in fig. 9. 
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Figure 4: Comparison between POWHEG and MCSNLO for the rapidity, invariant mass and transverse- 
momentum distributions of a Higgs boson with nin — 120 GeV, at the LHC pp coUider. 




Figure 5: Comparison between POWHEG and MCONLO for the rapidity of the leading jet and the ra- 
pidity difference of the Higgs boson and the leading jet, defined according to the SISCONE algorithm, 
with different jet cuts. 



4.2 POWHEG - PYTHIA comparison 

We now compare POWHEG and PYTHIA. The Higgs boson production implementation in 
PYTHIA includes matrix-element corrections, so that the p-r distribution of the Higgs boson 
is accurate at large px- In our comparisons, we always normalize the PYTHIA results to the 
full NLO cross section of POWHEG. We use the new px-ordered shower defined in the PYEVNW 
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Figure 6: Comparison between POWHEG and MCONLD for the transverse-momentum distributions of 
the leading jet, defined according to the SISCONE and the KT algorithms. 
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Figure 7: Comparison between POWHEG and MCONLO for the rapidity, invariant mass and transverse- 
momentum distributions of a Higgs boson with mn = 400 GeV, at the LHC pp collider. 



routine of PYTHIA, that should be more appropriate when interfacing to POWHEG. 

The only difference with respect to the POWHEG-MCONLO comparisons is in the generation 
of the Higgs boson virtuality, distributed now according to 

^ (M2 -mjjf + (M2 FH/mHy ' 
which is very similar to the form used in PYTHIA, except for the fact that PYTHIA includes 
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Figure 8: Comparison between POWHEG and MCONLO for the rapidity of the leading jet and the ra- 
pidity difference of the Higgs boson and the leading jet, defined according to the SISCONE algorithm, 
with different jet cuts. 
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Figure 9: Comparison of POWHEG, MCONLO and HERWIG (without matrix- element corrections), for 
the rapidity of the leading jet and the rapidity difference of the Higgs boson and the leading jet, 
defined according to the SISCONE algorithm, with different jet cuts. 



threshold effects in the calculation of the Higgs boson width. In fact, PYTHIA uses a running 
THiM^), that increases when a decay channel opens up. The effects of using a fixed or a 
running Th are more evident for a heavy Higgs boson, as will be shown in the following. 

In figs. 10 through 12 we compare results for the Tevatron pp collider, while in figs. 13 
through 15 we present results for the LHC. In all the plots we have set mn = 120 GeV. 
Results are in an impressive good agreement, both for inclusive quantities and for more 
exclusive ones. The only visible difference is in the transverse Higgs boson momentum dis- 
tribution at low at the LHC. This could be due to the different choice of the renormal- 
ization and factorization scale in the generation of radiation, our choice being constrained 
by the requirement of next-to-leading logarithmic accuracy in the Sudakov form factor. 

In fig. 16 we present a comparison with uih = 400 GeV. Mass thresholds effects in Th 
are evident in the invariant-mass distribution generated by PYTHIA. Below 2mt the total 
width is smaller than the fixed one we are using, and PYTHIA results are accordingly lower 
than ours. All other plots show instead good agreement with POWHEG. 
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Figure 10: Comparison between POWHEG and PYTHIA for the rapidity, invariant mass and 
transverse- momentum distributions of a Higgs boson with mu = 120 GeV, at Tevatron pp col- 
hder. PYTHIA outputs normahzed to the POWHEG cross section. 
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Figure 11: Comparison between POWHEG and PYTHIA for the rapidity of the leading jet and the ra- 
pidity difference of the Higgs boson and the leading jet, defined according to the SISCONE algorithm, 
with different jet cuts. PYTHIA outputs are normahzed to the POWHEG cross section. 



The good agreement between POWHEG and PYTHIA was to some extent expected. As 
already observed in refs. [9, 14], the matrix-element correction method used in PYTHIA [27, 
28] bears considerable similarities to POWHEG. 
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Figure 12: Comparison between POWHEG and PYTHIA for the transverse-momentum distributions 
of the leading jet, defined according to the SISCONE and the KT algorithms. PYTHIA outputs are 
normalized to the POWHEG cross section. 
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Figure 13: Comparison between POWHEG and PYTHIA for the rapidity, invariant mass and 
transverse-momentum distributions of a Higgs boson with mn = 120 GeV, at the LHC. PYTHIA 
outputs are normalized to the POWHEG cross section. 

4.3 The Pt distribution in POWHEG 

In this section we address the discrepancy in the px distributions in POWHEG and in MCONLO. 
First of all, we show in fig. 17 a comparison between the px spectrum of POWHEG, MCONLO 
and the NLO calculation. For sake of comparison, we have used in POWHEG and in the NLO 
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Figure 14: Comparison between POWHEG and PYTHIA for the rapidity of the leading jet and the ra- 
pidity difference of the Higgs boson and the leading jet, defined according to the SISCONE algorithm, 
with different jet cuts. PYTHIA outputs are normalized to the PDWHEG cross section. 




so 100 150 200 250 300 SO 100 150 200 250 300 

pi"* [GeV] pf [GeV] 

Figure 15: Comparison between POWHEG and PYTHIA for the transverse-momentum distributions 
of the leading jet, defined according to the SISCONE and the KT algorithms. PYTHIA outputs are 
normalized to the POWHEG cross section. 

calculation the same scale choice adopted in MCSNLO. We point out, however, that using a 
scale that depends upon the transverse momentum of radiation in POWHEG can only affect 
the B function. More specifically, one ends up using a transverse momentum dependent 
scale only in calculation of the real contributions in B, since the transverse momentum is 
zero for the Born, virtual and collinear remnant terms. Thus, this scale does not depend 
upon the transverse momentum of the real radiation, that is generated afterwards using 
the POWHEG Sudakov form factor. The choice of scale for radiation affects instead a single 
power of the coupling constant, since the Sudakov exponent is proportional to as- At low 
transverse momentum, this scale cannot be changed without spoiling the NLL accuracy of 
the Sudakov form factor. It can be changed, however, at large transverse momentum to 
explore further uncertainties. However, we have preferred not to implement this possibility. 
One should recall, in fact, that this scale only affects a single power of Og, and it thus has 
a much smaller effect than a scale change in the NLO cross section. 

We see from fig. 17 that MCONLO agrees better than POWHEG with the NLO calculation 
at large px- Since the difference between MCQNLO and POWHEG should be of next-to-next-to- 
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Figure 16: Comparison between POWHEG and PYTHIA for the rapidity, invariant mass and 
transverse- momentum distributions of a Higgs boson with niH — 400 GeV, at the LHC. PYTHIA 
outputs are normalized to the POWHEG cross section. 



leading order (NNLO), the difference between POWHEG and the NLO result should also be 
of NNLO. In fact we can easily trace the origin of this difference. From eq. (3.37), we infer 
that, at large the POWHEG differential cross section can be written as 



da 



(4.3) 



since the Sudakov form factor approaches 1 in this region. Neglecting the subdominant 
qq real contribution, this differs from the pure NLO result because of the presence of the 
factor 



1 + 0(03) 



(4.4) 



5(*i) 

It is known that radiative corrections in Higgs boson production are large, so that the 
0(03) term is in fact of order 1, and thus we find an enhancement that approaches a factor 
of two.'^ We have performed a clear cut test of this interpretation of the discrepancy. We 
have replaced the B function with the Born term B in the POWHEG program. The result of 



^We recall that normally the numerator and denominator in this factor are evaluated at different scales, 
since in B one uses a scale of the order of the Ifiggs boson transverse mass, while in the B term, one uses 
the transverse momentum. However, at large pt, the two scales become of the same order. 
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Figure 17: Comparison between POWHEG, MCONLO and tlic NLO calculation, for ttih = 120 GeV at 
the LHC. All calculations are performed in the mt — > oo approximation. Shower and hadronization 
are included in the MC results. The POWHEG result is also presented without shower and hadroniza- 
tion, and with a fixed-scale choice. 



this calculation is shown in comparison with the NLO curve in fig. 18. Since, as shown in 
fig. 17, the shower and hadronization are irrelevant for this distribution, we do not include 
them in the figure. In fig. 18 we have chosen to use independent renormalization and 
factorization scales, in order to perform a consistent comparison. Notice that, with this 
choice of scales, the NLO distribution is harder than the one shown in fig. 17. This is 
easily explained by the fact that the NLO process is proportional to (/in), and thus a p-r 
dependent renormalization scale can alter significantly the p^ distribution. 

At this point, we can ask whether the higher order terms included in POWHEG with the 
mechanism illustrated above do in fact give a reasonable estimate of true NNLO effects. 
We thus include in fig. 18 the NNLO result, obtained from the HMMLO program of ref. [16]. 
The result shows a rather good agreement between the NNLO result and POWHEG. Thus, 
our seemingly large corrections to the Higgs boson distributions are in fact very similar 
in size to the full NNLO result. Observe that in fig. 18 we have used a fixed scale choice 
for all the results. We were forced to do this, since the HNNLO program does not allow for 
other choices. However, because of the good agreement of the two POWHEG results in fig. 17, 
and because of the smaller scale dependence of the NNLO result, this should not make a 
severe difference. 

Because of a fortuitous circumstance, we did not need to worry about correcting for 
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Figure 18: Comparison between POWHEG and fixed NLO and NNLO distributions for the transverse- 
momentum of the Higgs boson. Plots are done for ttih = 120 GeV at the LHC. 



the large difference between the POWHEG and the NLO result at large radiation transverse 
momentum, since the known NNLO result seems to support the POWHEG one. We remark, 
however, that, had this not been the case, it is very easy to modify the POWHEG algorithm 
so to obtain a spectrum that agrees with the NLO calculation at large p^- This can be 
done as follows. Instead of using the full real cross section for the computation of the B 
function and of the Sudakov form factor, we can instead use a reduced real contribution 

R'^'^ = RxF, (4.5) 

where is a function of the real phase space, with F < 1 everywhere, such that F 
approaches 1 for small transverse momenta, and approaches zero for large transverse mo- 
menta. We perform the POWHEG generation using R^^'^ instead of R, and treat the remaining 
R X (1 — F) contribution to the cross section with the same method that we used for the 
Rqq contribution. This can be done, since i? x (1 — F) is dumped by the 1 — F factor in 
the singular region. It will then follow that, for large transverse momentum, the result 
would agree with the NLO calculation, since it would be dominated by the R x (1 — F) 
contribution. It turned out that, in all previous implementations, it was not necessary to 
use such procedure. As remarked before, thanks to the known properties of the NNLO 
result, this was not necessary even in this case. We have however performed such study, 
just in order to illustrate the flexibility of the POWHEG method. We have chosen for F the 
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following form 



F 



(4.6) 



The resulting transverse-momentum distribution at the LHC, for a Higgs boson mass of 
400 GeV, is shown in fig. 19 for /i ^ oo (standard POWHEG), h = 120 GeV and h = 400 GeV. 
One can see that it is not difficult to get distributions that undershoot the MCONLO one in 
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Figure 19: Comparison of the predictions of MCONLO, standard POWHEG {h — > oo) and POWHEG with 
two different values of the parameter h {h = 120 GeV and h = uriH — 400 GeV) in the function F 
of eq. (4.6), for the transverse-momentum distributions of a Higgs boson, at the LHC pp collider. 



the intermediate range of px- We also observe that, with this procedure, no undesired 
features of other distributions appear. In particular, the distribution in the rapidity of the 
hardest jet, and in the rapidity difference between the hardest jet and the Higgs boson 
remain qualitatively the same, as shown in fig. 20. 

4.4 Next-to-leading logarithmic resummation 

As explained in section (4.4) of ref. [9], one can reach next-to-leading logarithmic (NLL) 
accuracy of soft gluon resummation if the number of coloured partons involved in the hard 
scattering is less or equal to three. This can be obtained by replacing the strong coupling 
constant in the Sudakov exponent with [29] 
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Figure 20: Comparison of the predictions of the standard POWHEG {h oo), and PDWHEG with two 
different values of the parameter h [h ^ 120 GeV and h = ran = 400 GcV) in the function F of 
eq. (4.6), for the rapidity of the leading jet and the rapidity difference of the Higgs boson and the 
leading jet, defined according to the SISCONE algorithm, with different jet cuts, at the LHC. 
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Figure 21: Comparison between POWHEG and HqT for the transverse-momentum distributions of a 
Higgs boson, at the LHC. The POWHEG result without shower and hadronization is also shown. 



where the MS, 1-loop expression of as should be used. The previous replacement may 
also be implemented by a simple redefinition of the strong scale A, which, for five active 
flavours (nf = 5), becomes Amc = 1.569 A^. We have exploited this possibility in our 
code, so that our result should agree with the NLL resummed one. A comparison has 
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been thus carried out with the HqT [30] program, that performs such a resummation. We 
have adopted fixed renormahzation and factorization scales. Results arc shown in fig. 21, 
together with next-to-next-to-leading logarithmic (NNLL) resummation, always from HqT, 
just for reference purposes. We see a fair agreement between the POWHEG result and the NLL 
analytic one, as expected. The different behaviour of the POWHEG result without shower 
and hadronization at very low may be ascribed to the particular implementation of the 
minimum transverse momentum that we use, that is, to a large extent, arbitrary. 

We observe that, in all cases, we do not expect full agreement between the POWHEG 
result without shower, and the NLL calculation. In fact, the POWHEG curve without shower 
represents the Sudakov form factor for the Pt of the hardest emission, while, in the NLL 
calculation, the total p-^ distribution (i.e. the sum of the transverse momenta of all emis- 
sions) is considered. Thus, it is only after the inclusion of the full shower effects that the 
two distributions have a meaningful comparison. 

5. Conclusions 

In this paper we have reported on a complete implementation of Higgs boson production via 
gluon fusion at next-to-leading order in QCD, in the POWHEG framework. The calculation 
was performed within the Prixione-Kunszt-Signer [20, 21] subtraction approach. We have 
also shown how to deal with non-singular real contributions, that do not present a valid 
underlying Born matrix element. 

The results of our work have been compared extensively with MCONLO and PYTHIA 
shower Monte Carlo programs. The PYTHIA results, normalized to the total NLO cross 
section, are in good agreement with POWHEG, except for differences in the low transverse- 
momentum distributions of the Higgs boson at the LHC. The MCSNLO results are in fair 
agreement with POWHEG, except for the px distribution of the Higgs boson, and consequently 
of the hardest jet, in the high-px region. In this region the POWHEG distributions are 
generally harder. We have shown that this is due to NNLO effects in the POWHEG formula 
for the differential cross section. We checked that these effects actually bring our result 
closer to the NNLO one [16]. The low-px region was instead tested against the analytic 
resummed results [30]. We find again good agreement up to NLL accuracy. 

Furthermore, we have also examined the distributions in the difference of the hardest 
jet and the Higgs boson rapidity. The dip found in previous implementations [10, 11, 14] 
is still present. We remark that this seems to be a general feature of MCONLO, since other 
calculations do not find effects of this kind [17, 18, 19]. 

The computer code for the POWHEG implementations presented in this paper is avail- 
able, together with the manual, at the site 

http : //moby .mib . inf n . it/~nason/POWHEG 
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